sc-RNAseq analysis of breast cancer
Alba Machado-Lopez
October 2022
Load libraries
if (!require(R.utils)) install.packages('R.utils')## Warning: package 'R.utils' was built under R version 4.0.5library(R.utils)
if (!require(dplyr)) install.packages('dplyr')## Warning: package 'dplyr' was built under R version 4.0.5library(dplyr)
if (!require(stringr)) install.packages('stringr')
library(stringr)
if (!require(Seurat)) install.packages('Seurat')
library(Seurat)
if (!require(purrr)) install.packages('purrr')
library(purrr)
if (!require(scCATCH)) install.packages('scCATCH')
library(scCATCH)
if (!require(ggplot2)) install.packages('ggplot2')## Warning: package 'ggplot2' was built under R version 4.0.5library(ggplot2)
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!require(SingleR)) install.packages('SingleR')## Warning: package 'matrixStats' was built under R version 4.0.5## Warning: package 'BiocGenerics' was built under R version 4.0.5## Warning: package 'GenomeInfoDb' was built under R version 4.0.5library(SingleR)
if (!require(msigdbr)) install.packages('msigdbr')## Warning: package 'msigdbr' was built under R version 4.0.5library(msigdbr)
if (!require(tibble)) install.packages('tibble')## Warning: package 'tibble' was built under R version 4.0.5library(tibble)
if (!require(fgsea)) install.packages('fgsea')
library(fgsea)
#if (!require(clusterProfiler)) install.packages('clusterProfiler')
library(clusterProfiler)
#if (!require(org.Hs.eg.db)) install.packages('org.Hs.eg.db')
library(org.Hs.eg.db)
if (!require("org.Mm.eg.db", quietly = TRUE))
BiocManager::install("org.Mm.eg.db")
library(org.Mm.eg.db)Obtain files
Download data
This project uses the data from Grosselin et al, 2019, which is uploaded to GEO.
mkdir -p data
cd data
wget -r -np -nd 'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE117nnn/GSE117309/suppl/GSE117309_RAW.tar' -R "index.html*" # Omitting the index fileDecompress files
# Untar the main file
untar("data/GSE117309_RAW.tar", exdir="data/")
# Untar sc-RNAseq files
RNA_files <- list.files(path="data/",pattern="scRNA")
untar_RNA_files <- function(file){untar(paste0("data/",file), exdir="data/RNA_files")}
sapply(RNA_files, untar_RNA_files)Read files into R
# Create function to load the files
load_scRNA_files <- function(model, organism){
directory <- paste0("data/RNA_files/filtered_gene_bc_matrices_",model,"/",organism,"/")
expression_matrix <- ReadMtx(
mtx=paste0(directory, "matrix.mtx"),
features=paste0(directory, "genes.tsv"),
cells=paste0(directory, "barcodes.tsv"))
# Correct gene names for easier automatization
expression_matrix@Dimnames[[1]] <- str_remove_all(expression_matrix@Dimnames[[1]],
pattern=paste0(organism,"_"))
seurat_object <- CreateSeuratObject(counts = expression_matrix, project=model)
return(seurat_object)
}
models <- c("HBCx-95","HBCx-95_CAPAR", "HBCx-22","HBCx22-TAMR")
scRNA_files_hg19 <- mapply(FUN=load_scRNA_files, models, organism="hg19")## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')names(scRNA_files_hg19) <- paste0(names(scRNA_files_hg19), "_hg19")
scRNA_files_mm10 <- mapply(FUN=load_scRNA_files, models, organism="mm10")
names(scRNA_files_mm10) <- paste0(names(scRNA_files_mm10), "_mm10")Analysis of human cells
Data integration
# Normalize and identify variable features
scRNA_files_hg19 <- lapply(X = scRNA_files_hg19, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
# Select features repeatedly variable between the 4 datasets
features <- SelectIntegrationFeatures(object.list = scRNA_files_hg19)
# Select anchors for integration
brca_anchors <- FindIntegrationAnchors(object.list = scRNA_files_hg19, anchor.features = features)## Warning in CheckDuplicateCellNames(object.list = object.list): Some cell names
## are duplicated across objects provided. Renaming to enforce unique cell names.## Scaling features for provided objects## Finding all pairwise anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 2276 anchors## Filtering anchors## Retained 1594 anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 2558 anchors## Filtering anchors## Retained 1046 anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 3666 anchors## Filtering anchors## Retained 1570 anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 2201 anchors## Filtering anchors## Retained 1215 anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 2766 anchors## Filtering anchors## Retained 1635 anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 3104 anchors## Filtering anchors## Retained 2412 anchors# Create the integrated data assay
scRNA_combined <- IntegrateData(anchorset=brca_anchors)## Merging dataset 4 into 3## Extracting anchors for merged samples## Finding integration vectors## Finding integration vector weights## Integrating data## Merging dataset 1 into 2## Extracting anchors for merged samples## Finding integration vectors## Finding integration vector weights## Integrating data## Merging dataset 2 1 into 3 4## Extracting anchors for merged samples## Finding integration vectors## Finding integration vector weights## Integrating dataQuality control
# Use the integrated assay as default
DefaultAssay(scRNA_combined) <- "RNA"
# Creating a new slot with the percentage of counts on mitochondrial reads
scRNA_combined[["percent.mt"]] <- PercentageFeatureSet(scRNA_combined, pattern="MT-")
# Violin plots of specific features
VlnPlot(scRNA_combined,
features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
ncol = 3,
split.by="orig.ident")## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##
## This message will be shown once per session.
We will remove cells with a percentage of mitochondrial reads too high (>5%) and with an abnormally high or low number of feature counts, that can suggest empty droplets or multiplets in the microfluidic system, respectively.
scRNA_combined <- subset(scRNA_combined, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 10)
VlnPlot(scRNA_combined,
features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
ncol = 3,
split.by="orig.ident")
Visualization and clustering
Check integration
# Change the default to the integrated dataset
DefaultAssay(scRNA_combined) <- "integrated"
scRNA_combined <- ScaleData(scRNA_combined, verbose = FALSE)
scRNA_combined <- RunPCA(scRNA_combined, npcs = 30, verbose = FALSE)
scRNA_combined <- RunUMAP(scRNA_combined, reduction = "pca", dims = 1:30)## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 22:19:16 UMAP embedding parameters a = 0.9922 b = 1.112## 22:19:16 Read 3006 rows and found 30 numeric columns## 22:19:16 Using Annoy for neighbor search, n_neighbors = 30## 22:19:16 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 22:19:17 Writing NN index file to temp file C:\Users\alba_\AppData\Local\Temp\RtmpUB6mHX\file2c705eea473b
## 22:19:17 Searching Annoy index using 1 thread, search_k = 3000
## 22:19:18 Annoy recall = 100%
## 22:19:19 Commencing smooth kNN distance calibration using 1 thread
## 22:19:21 Initializing from normalized Laplacian + noise
## 22:19:22 Commencing optimization for 500 epochs, with 129862 positive edges
## 22:19:34 Optimization finishedscRNA_combined <- FindNeighbors(scRNA_combined, reduction = "pca", dims = 1:30)## Computing nearest neighbor graph
## Computing SNNscRNA_combined <- FindClusters(scRNA_combined, resolution = 0.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 3006
## Number of edges: 151626
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7945
## Number of communities: 5
## Elapsed time: 0 seconds# Visualization
p1 <- DimPlot(scRNA_combined, reduction = "umap", group.by = "orig.ident")
p2 <- DimPlot(scRNA_combined, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2
DimPlot(scRNA_combined, reduction = "umap", split.by = "orig.ident")
After normalization and integration, the UMAPs of all subsets look quite similar, so we can assume integration is correct and we can proceed forward. Once this has been verified, we can go back a few steps and refine our analysis
Determine the dimensionality of the dataset
scRNA_combined <- JackStraw(scRNA_combined, num.replicate = 100)
scRNA_combined <- ScoreJackStraw(scRNA_combined, dims=1:20)
JackStrawPlot(scRNA_combined, dims=1:20)## Warning: Removed 28000 rows containing missing values (geom_point).
ElbowPlot(scRNA_combined)
Based on these two graphs, we may consider using 13-15 PCs.
scRNA_combined <- FindNeighbors(scRNA_combined, reduction = "pca", dims = 1:13)## Computing nearest neighbor graph## Computing SNNscRNA_combined <- FindClusters(scRNA_combined, resolution = 0.6) # Try different levels of resolution and chose the one that results in cluster with most biological sense1## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 3006
## Number of edges: 112761
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7863
## Number of communities: 6
## Elapsed time: 0 secondsscRNA_combined <- RunUMAP(scRNA_combined, dims = 1:13)## 22:22:57 UMAP embedding parameters a = 0.9922 b = 1.112## 22:22:57 Read 3006 rows and found 13 numeric columns## 22:22:57 Using Annoy for neighbor search, n_neighbors = 30## 22:22:57 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 22:22:58 Writing NN index file to temp file C:\Users\alba_\AppData\Local\Temp\RtmpUB6mHX\file2c7055085643
## 22:22:58 Searching Annoy index using 1 thread, search_k = 3000
## 22:23:00 Annoy recall = 100%
## 22:23:01 Commencing smooth kNN distance calibration using 1 thread
## 22:23:03 Initializing from normalized Laplacian + noise
## 22:23:03 Commencing optimization for 500 epochs, with 123068 positive edges
## 22:23:17 Optimization finishedDimPlot(scRNA_combined, reduction="umap")
DimPlot(scRNA_combined, reduction = "umap", split.by = "orig.ident")
Finding marker genes
hg19_markers <- FindAllMarkers(scRNA_combined, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)## Calculating cluster 0## Calculating cluster 1## Calculating cluster 2## Calculating cluster 3## Calculating cluster 4## Calculating cluster 5hg19_markers %>%
group_by(cluster) %>%
slice_max(n = 2, order_by = avg_log2FC)## # A tibble: 12 x 7
## # Groups: cluster [6]
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
## 1 3.69e- 69 0.815 0.926 0.82 7.38e- 66 0 CLEC3A
## 2 3.34e- 72 0.767 0.939 0.843 6.69e- 69 0 AGR2
## 3 1.15e- 79 0.879 0.997 0.931 2.30e- 76 1 NDRG1
## 4 3.76e- 30 0.500 0.985 0.889 7.51e- 27 1 IGFBP5
## 5 2.28e- 6 1.55 0.716 0.73 4.57e- 3 2 SCGB2A2
## 6 3.31e- 88 1.11 0.973 0.896 6.62e- 85 2 IGFBP5
## 7 1.53e-161 2.03 0.969 0.7 3.06e-158 3 HIST1H4C
## 8 2.69e-218 1.93 0.971 0.632 5.38e-215 3 KIAA0101
## 9 5.22e-142 3.11 0.951 0.561 1.04e-138 4 KRT16
## 10 5.35e-135 2.87 0.904 0.447 1.07e-131 4 KRT6B
## 11 6.67e- 90 2.70 0.98 0.71 1.33e- 86 5 UBE2C
## 12 8.13e-113 2.28 0.995 0.686 1.63e-109 5 PTTG1hg19_markers %>%
group_by(cluster) %>%
top_n(n = 10, wt = avg_log2FC) -> top10
DoHeatmap(scRNA_combined, features = top10$gene) + NoLegend()
Functional analysis
Gene set enrichment analysis
MSigDB_df <- msigdbr(species = "Homo sapiens")
fgsea_sets<- MSigDB_df %>% split(x = .$gene_symbol, f = .$gs_name)
plot_GSEA <- function(cluster_number){
genes <- hg19_markers %>% filter(cluster==as.character(cluster_number)) %>%
arrange(desc(p_val_adj)) %>%
dplyr::select(gene, avg_log2FC)
# Create dataframe with genes
ranks <- deframe(genes)
# Perform GSEA analysis
fgseaRes <- fgsea(fgsea_sets,
stats = ranks,
minSize=10,
maxSize=500,
nperm=1000000)
# Plot
ggplot(fgseaRes %>% filter(padj < 0.01) %>% head(n= 20), aes(reorder(pathway, NES), NES)) +
geom_col(aes(fill= NES < 7.5)) +
coord_flip() +
labs(x="Pathway",
y="Normalized Enrichment Score",
title=paste0("Hallmark pathways NES from GSEA for cluster",as.character(cluster_number)))
}
for (cluster in 0:5){plot_GSEA(cluster)}Gene ontology
plot_GO <- function(cluster_number){
genes <- hg19_markers %>% filter(cluster==as.character(cluster_number))
# Calculate enrichment
enrichment <- enrichGO(gene = genes$gene,
OrgDb = org.Hs.eg.db,
keyType = "SYMBOL",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
enrichment <- clusterProfiler::simplify(enrichment, cutoff=0.7, by="p.adjust", select_fun=min)
GO_ggdata <- enrichment %>%
as_data_frame() %>%
arrange(Count)
GO_ggdata$Description <- factor(GO_ggdata$Description, levels = GO_ggdata$Description)
ggplot(GO_ggdata, aes(x = Description, y = Count, fill = p.adjust)) +
geom_bar(stat = "identity") +
scale_colour_viridis_d(begin=0,end=1) +
coord_flip() +
ylab("Number of genes") +
xlab("GO Terms") +
theme(axis.text.y = element_text(size=10))
}
for (cluster in 0:5){print(plot_GO(cluster))}## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
Automated cluster annotation
Method 1- scCATCH
hg19_geneinfo <- rev_gene(data = scRNA_combined[['integrated']]@data,
data_type = "data",
species="Human",
geneinfo = geneinfo)
hg19_ann_1 <- createscCATCH(data = scRNA_combined[['integrated']]@data,
cluster = as.character(Idents(scRNA_combined)))
hg19_ann_1 <- findmarkergene(object = hg19_ann_1, species = "Human", marker = cellmatch, tissue = "Breast", cancer="Breast Cancer", use_method = "2")## There are 535 potential marker genes in CellMatch database for Human Breast Cancer on Breast.hg19_ann_1 <- findcelltype(hg19_ann_1)
hg19_ann_1@celltype## cluster cluster_marker
## 1 1 FN1, IGFBP5, TACSTD2, PMEPA1, CD44
## 2 2 IGFBP5
## 3 3 CENPF, BIRC5, HELLS, BLVRB
## 4 0 BST2
## 5 5 CENPF, BIRC5
## 6 4 G0S2, FN1, ITGA6, TACSTD2, LIMA1, PMEPA1, DNAJB1, ZFP36L2, CD44
## cell_type celltype_score celltype_related_marker
## 1 Cancer Stem Cell 0.98 CD44
## 2 Killer Cell 0.50 IGFBP5
## 3 Helper T Cell 0.61 CENPF, BIRC5, HELLS
## 4 T Cell 0.50 BST2
## 5 Helper T Cell 0.58 CENPF, BIRC5
## 6 Cancer Stem Cell 0.98 CD44, ITGA6
## PMID
## 1 24099815, 28986882, 27840965, 23799994, 23053657, 23768049, 28719381, 28579829, 28554844, 28399903, 28166194, 28148288, 28036386, 28012234, 27768764, 27630305, 27458252, 27133070, 27115469, 26893363, 26170011, 25990436, 24596390, 24297508, 24171482, 24144739, 24055090, 23543868, 23112116, 23056395, 23036368, 22761812, 22464177, 21979823, 21802218, 21092249, 22459349, 27640754, 25940879
## 2 28474673
## 3 28474673
## 4 28474673
## 5 28474673
## 6 24099815, 28986882, 27840965, 23799994, 23053657, 23768049, 28719381, 28579829, 28554844, 28399903, 28166194, 28148288, 28036386, 28012234, 27768764, 27630305, 27458252, 27133070, 27115469, 26893363, 26170011, 25990436, 24596390, 24297508, 24171482, 24144739, 24055090, 23543868, 23112116, 23056395, 23036368, 22761812, 22464177, 21979823, 21802218, 21092249, 22459349, 27640754, 25940879Method 2- SingleR
hpca.se <- HumanPrimaryCellAtlasData()## Warning: 'HumanPrimaryCellAtlasData' is deprecated.
## Use 'celldex::HumanPrimaryCellAtlasData' instead.
## See help("Deprecated")## snapshotDate(): 2020-10-27## see ?celldex and browseVignettes('celldex') for documentation## loading from cache## see ?celldex and browseVignettes('celldex') for documentation## loading from cachehpca.se## class: SummarizedExperiment
## dim: 19363 713
## metadata(0):
## assays(1): logcounts
## rownames(19363): A1BG A1BG-AS1 ... ZZEF1 ZZZ3
## rowData names(0):
## colnames(713): GSM112490 GSM112491 ... GSM92233 GSM92234
## colData names(3): label.main label.fine label.onthg19_ann_2 <- SingleR(test = scRNA_combined[['integrated']]@data,
ref = hpca.se,
assay.type.test=1,
labels = hpca.se$label.main)
hg19_ann_2## DataFrame with 3006 rows and 5 columns
## scores first.labels
## <matrix> <character>
## AAACCTGGTTCCACTC-1_1 0.284388:0.263078:0.208514:... Epithelial_cells
## AACCGCGCACCTTGTC-1_1 0.219088:0.219750:0.185685:... Epithelial_cells
## AACCGCGTCCAAGCCG-1_1 0.288097:0.332301:0.294157:... Endothelial_cells
## AACGTTGGTTCAGTAC-1_1 0.312886:0.359925:0.326496:... BM & Prog.
## AACGTTGTCCTGCTTG-1_1 0.254901:0.215151:0.166806:... Epithelial_cells
## ... ... ...
## TTCTTAGCATTCTTAC-1_4 0.128125:0.109205:0.078661:... Hepatocytes
## TTGACTTGTATGCTTG-1_4 0.171581:0.194234:0.169912:... Hepatocytes
## TTGACTTTCCCACTTG-1_4 0.175745:0.202312:0.173485:... Epithelial_cells
## TTGGCAAGTCTTGTCC-1_4 0.192326:0.156040:0.144870:... Epithelial_cells
## TTTACTGGTAGGCATG-1_4 0.140242:0.201939:0.168875:... Epithelial_cells
## tuning.scores labels pruned.labels
## <DataFrame> <character> <character>
## AAACCTGGTTCCACTC-1_1 0.223729: 0.2155909 Epithelial_cells Epithelial_cells
## AACCGCGCACCTTGTC-1_1 0.244705:-0.0810895 Epithelial_cells Epithelial_cells
## AACCGCGTCCAAGCCG-1_1 0.183081: 0.1802609 Endothelial_cells Endothelial_cells
## AACGTTGGTTCAGTAC-1_1 0.297462: 0.2850939 BM & Prog. BM & Prog.
## AACGTTGTCCTGCTTG-1_1 0.220434: 0.1080830 Epithelial_cells Epithelial_cells
## ... ... ... ...
## TTCTTAGCATTCTTAC-1_4 0.1070612:0.0436891 Epithelial_cells Epithelial_cells
## TTGACTTGTATGCTTG-1_4 0.0747311:0.0635165 Hepatocytes Hepatocytes
## TTGACTTTCCCACTTG-1_4 0.3865267:0.3086367 Epithelial_cells Epithelial_cells
## TTGGCAAGTCTTGTCC-1_4 0.2565691:0.1856667 Epithelial_cells Epithelial_cells
## TTTACTGGTAGGCATG-1_4 0.1900828:0.1195525 Epithelial_cells Epithelial_cellsAnalysis of mouse cells
Data integration
# Normalize and identify variable features
scRNA_files_mm10 <- lapply(X = scRNA_files_mm10, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
# Select features repeatedly variable between the 4 datasets
features_mm10 <- SelectIntegrationFeatures(object.list = scRNA_files_mm10)
# Select anchors for integration
brca_anchors_mm10 <- FindIntegrationAnchors(object.list = scRNA_files_mm10, anchor.features = features_mm10)## Warning in CheckDuplicateCellNames(object.list = object.list): Some cell names
## are duplicated across objects provided. Renaming to enforce unique cell names.## Scaling features for provided objects## Finding all pairwise anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 2294 anchors## Filtering anchors## Retained 2030 anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 2088 anchors## Filtering anchors## Retained 1826 anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 2833 anchors## Filtering anchors## Retained 2294 anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 2231 anchors## Filtering anchors## Retained 2108 anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 2983 anchors## Filtering anchors## Retained 2716 anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 2769 anchors## Filtering anchors## Retained 2515 anchors# Create the integrated data assay
scRNA_combined_mm10 <- IntegrateData(anchorset=brca_anchors_mm10)## Merging dataset 1 into 4## Extracting anchors for merged samples## Finding integration vectors## Finding integration vector weights## Integrating data## Merging dataset 2 into 4 1## Extracting anchors for merged samples## Finding integration vectors## Finding integration vector weights## Integrating data## Merging dataset 3 into 4 1 2## Extracting anchors for merged samples## Finding integration vectors## Finding integration vector weights## Integrating dataQuality control
# Use the integrated assay as default
DefaultAssay(scRNA_combined_mm10) <- "RNA"
# Creating a new slot with the percentage of counts on mitochondrial reads
scRNA_combined_mm10[["percent.mt"]] <- PercentageFeatureSet(scRNA_combined_mm10, pattern="mt-")
# Violin plots of specific features
VlnPlot(scRNA_combined_mm10,
features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
ncol = 3,
split.by="orig.ident")
We will remove cells with a percentage of mitochondrial reads too high (>5%) and with an abnormally high or low number of feature counts, that can suggest empty droplets or multiplets in the microfluidic system, respectively.
scRNA_combined_mm10 <- subset(scRNA_combined_mm10, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 10)
VlnPlot(scRNA_combined_mm10,
features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
ncol = 3,
split.by="orig.ident")
Visualization and clustering
Check integration
# Change the default to the integrated dataset
DefaultAssay(scRNA_combined_mm10) <- "integrated"
scRNA_combined_mm10 <- ScaleData(scRNA_combined_mm10, verbose = FALSE)
scRNA_combined_mm10 <- RunPCA(scRNA_combined_mm10, npcs = 30, verbose = FALSE)
scRNA_combined_mm10 <- RunUMAP(scRNA_combined_mm10, reduction = "pca", dims = 1:30)## 22:31:28 UMAP embedding parameters a = 0.9922 b = 1.112## 22:31:28 Read 4017 rows and found 30 numeric columns## 22:31:28 Using Annoy for neighbor search, n_neighbors = 30## 22:31:28 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 22:31:29 Writing NN index file to temp file C:\Users\alba_\AppData\Local\Temp\RtmpUB6mHX\file2c70633baa4
## 22:31:29 Searching Annoy index using 1 thread, search_k = 3000
## 22:31:31 Annoy recall = 100%
## 22:31:33 Commencing smooth kNN distance calibration using 1 thread
## 22:31:35 Initializing from normalized Laplacian + noise
## 22:31:36 Commencing optimization for 500 epochs, with 167364 positive edges
## 22:31:56 Optimization finishedscRNA_combined_mm10 <- FindNeighbors(scRNA_combined_mm10, reduction = "pca", dims = 1:30)## Computing nearest neighbor graph
## Computing SNNscRNA_combined_mm10 <- FindClusters(scRNA_combined_mm10, resolution = 0.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 4017
## Number of edges: 155514
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9021
## Number of communities: 14
## Elapsed time: 0 seconds# Visualization
p1 <- DimPlot(scRNA_combined_mm10, reduction = "umap", group.by = "orig.ident")
p2 <- DimPlot(scRNA_combined_mm10, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2
DimPlot(scRNA_combined_mm10, reduction = "umap", split.by = "orig.ident")
After normalization and integration, the UMAPs of all subsets look quite similar, so we can assume integration is correct and we can proceed forward. Once this has been verified, we can go back a few steps and refine our analysis
Determine the dimensionality of the dataset
scRNA_combined_mm10 <- JackStraw(scRNA_combined_mm10, num.replicate = 100)
scRNA_combined_mm10 <- ScoreJackStraw(scRNA_combined_mm10, dims=1:20)
JackStrawPlot(scRNA_combined_mm10, dims=1:20)## Warning: Removed 28000 rows containing missing values (geom_point).
ElbowPlot(scRNA_combined_mm10)
Based on these two graphs, we may consider using 10 PCs.
scRNA_combined_mm10 <- FindNeighbors(scRNA_combined_mm10, reduction = "pca", dims = 1:10)## Computing nearest neighbor graph## Computing SNNscRNA_combined_mm10 <- FindClusters(scRNA_combined_mm10, resolution = 0.6) # Try different levels of resolution and chose the one that results in cluster with most biological sense1## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 4017
## Number of edges: 128730
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8992
## Number of communities: 12
## Elapsed time: 0 secondsscRNA_combined_mm10 <- RunUMAP(scRNA_combined_mm10, dims = 1:10)## 22:36:02 UMAP embedding parameters a = 0.9922 b = 1.112## 22:36:02 Read 4017 rows and found 10 numeric columns## 22:36:02 Using Annoy for neighbor search, n_neighbors = 30## 22:36:02 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 22:36:03 Writing NN index file to temp file C:\Users\alba_\AppData\Local\Temp\RtmpUB6mHX\file2c7063c14705
## 22:36:03 Searching Annoy index using 1 thread, search_k = 3000
## 22:36:05 Annoy recall = 100%
## 22:36:07 Commencing smooth kNN distance calibration using 1 thread
## 22:36:10 Initializing from normalized Laplacian + noise
## 22:36:10 Commencing optimization for 500 epochs, with 160052 positive edges
## 22:36:29 Optimization finishedDimPlot(scRNA_combined_mm10, reduction="umap")
DimPlot(scRNA_combined_mm10, reduction = "umap", split.by = "orig.ident")
Finding marker genes
mm10_markers <- FindAllMarkers(scRNA_combined_mm10, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)## Calculating cluster 0## Calculating cluster 1## Calculating cluster 2## Calculating cluster 3## Calculating cluster 4## Calculating cluster 5## Calculating cluster 6## Calculating cluster 7## Calculating cluster 8## Calculating cluster 9## Calculating cluster 10## Calculating cluster 11mm10_markers %>%
group_by(cluster) %>%
slice_max(n = 2, order_by = avg_log2FC)## # A tibble: 24 x 7
## # Groups: cluster [12]
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
## 1 4.20e-194 2.98 0.868 0.692 8.39e-191 0 Gzme
## 2 3.12e-183 2.86 0.96 0.748 6.24e-180 0 Spp1
## 3 1.17e-304 3.83 0.994 0.742 2.35e-301 1 C1qa
## 4 0 3.61 1 0.725 0 1 Ctss
## 5 6.73e-127 2.38 0.871 0.597 1.35e-123 2 Mmp3
## 6 2.57e-223 2.10 0.945 0.547 5.14e-220 2 Mfap4
## 7 6.65e-190 2.73 0.976 0.647 1.33e-186 3 Apod
## 8 3.04e-184 2.55 0.988 0.797 6.08e-181 3 Gsn
## 9 6.90e-169 3.22 0.997 0.696 1.38e-165 4 Clec3b
## 10 1.34e-169 3.21 0.983 0.409 2.67e-166 4 Sbsn
## # ... with 14 more rowsmm10_markers %>%
group_by(cluster) %>%
top_n(n = 10, wt = avg_log2FC) -> top10
DoHeatmap(scRNA_combined_mm10, features = top10$gene) + NoLegend()
Functional analysis
Gene set enrichment analysis
MSigDB_df <- msigdbr(species = "Homo sapiens")
fgsea_sets<- MSigDB_df %>% split(x = .$gene_symbol, f = .$gs_name)
plot_GSEA <- function(cluster_number){
genes <- mm10_markers %>% filter(cluster==as.character(cluster_number)) %>%
arrange(desc(p_val_adj)) %>%
dplyr::select(gene, avg_log2FC)
# Create dataframe with genes
ranks <- deframe(genes)
# Perform GSEA analysis
fgseaRes <- fgsea(fgsea_sets,
stats = ranks,
minSize=10,
maxSize=500,
nperm=1000000)
# Plot
ggplot(fgseaRes %>% filter(padj < 0.01) %>% head(n= 20), aes(reorder(pathway, NES), NES)) +
geom_col(aes(fill= NES < 7.5)) +
coord_flip() +
labs(x="Pathway",
y="Normalized Enrichment Score",
title=paste0("Hallmark pathways NES from GSEA for cluster",as.character(cluster_number)))
}
plot_GSEA(0)
for (cluster in 0:5){plot_GSEA(cluster)}Gene ontology
plot_GO_mm10 <- function(cluster_number){
genes <- mm10_markers %>% filter(cluster==as.character(cluster_number))
# Calculate enrichment
enrichment <- enrichGO(gene = genes$gene,
OrgDb = org.Mm.eg.db,
keyType = "SYMBOL",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
enrichment <- clusterProfiler::simplify(enrichment, cutoff=0.7, by="p.adjust", select_fun=min)
GO_ggdata <- enrichment %>%
as_data_frame() %>%
arrange(Count)
GO_ggdata$Description <- factor(GO_ggdata$Description, levels = GO_ggdata$Description)
ggplot(GO_ggdata, aes(x = Description, y = Count, fill = p.adjust)) +
geom_bar(stat = "identity") +
scale_colour_viridis_d(begin=0,end=1) +
coord_flip() +
ylab("Number of genes") +
xlab("GO Terms") +
theme(axis.text.y = element_text(size=10))
}
for (cluster in 0:5){print(plot_GO_mm10(cluster))}
Automated cluster annotation
Method 1- scCATCH
MouseRNAseqData(ensembl = FALSE, cell.ont = c("all", "nonna", "none"))## snapshotDate(): 2020-10-27## see ?celldex and browseVignettes('celldex') for documentation## loading from cache## see ?celldex and browseVignettes('celldex') for documentation## loading from cache## class: SummarizedExperiment
## dim: 21214 358
## metadata(0):
## assays(1): logcounts
## rownames(21214): Xkr4 Rp1 ... LOC100039574 LOC100039753
## rowData names(0):
## colnames(358): ERR525589Aligned ERR525592Aligned ... SRR1044043Aligned
## SRR1044044Aligned
## colData names(3): label.main label.fine label.ontmm10_geneinfo <- rev_gene(data = scRNA_combined_mm10[['integrated']]@data,
data_type = "data",
species="Mouse",
geneinfo = geneinfo)
mm10_ann_1 <- createscCATCH(data = scRNA_combined_mm10[['integrated']]@data,
cluster = as.character(Idents(scRNA_combined_mm10)))
mm10_ann_1 <- findmarkergene(object = mm10_ann_1, species = "Mouse", marker = cellmatch, tissue="Mammary gland", use_method = "2")## There are 380 potential marker genes in CellMatch database for Mouse on Mammary gland.mm10_ann_1 <- findcelltype(mm10_ann_1)
mm10_ann_1@celltype## cluster
## 1 4
## 2 5
## 3 0
## 4 7
## 5 1
## 6 3
## 7 6
## 8 2
## 9 8
## 10 9
## 11 11
## 12 10
## cluster_marker
## 1 Igfbp5, Apod, Mt1, Clec3b, Gsn, Aldh1a3, Dcn, Plet1, Serpinb1a, Igfbp4, Efemp1, Procr, Igfbp6, Rarres2, Fn1, Htra3, S100a4, Ogn, Dpt, Ly6c1, Col3a1, Lum, Gadd45g, Ctsl, Mfap5, Serping1, Glul, Gpx3, Ly6a, Dpep1, Serpinf1, Ctsk, Pcolce, Loxl1, Fstl1, Phlda1, Plpp3, S100a6, Aebp1, Pdpn, C1s1, Mmp2, Lsp1, Anxa1, Anxa2
## 2 Fabp4, Ptn, Gng11, Lrg1, Col4a2, Atf3, Procr, Col4a1, Ly6c1, Id3, Rbp1, Fscn1, Crip2, Sparc, Csrp1, Tspan13, Ly6a, Cdkn1a, Tsc22d1, Kit, Itga6, Cd200, Gata2, Cnn3, Slc12a2, Igfbp7, S100a6, Serpinh1, Vim, Anxa2
## 3 Actg2, Gng11, Acta2, Igfbp2, Tagln, Ctsd, Myl9, Cnn1, Col4a2, Postn, Tpm2, Igfbp6, Tpm1, Col6a3, Col1a1, Col1a2, Fn1, Col3a1, Rbp1, Fscn1, Crip2, Col6a1, Mfap5, Sparc, Csrp1, Mfge8, Col6a2, Gpx3, Bgn, Lgals1, Tsc22d1, Serpinf1, Col5a2, Cd200, Cnn3, Pcolce, Loxl1, Fstl1, Mmp14, Igfbp7, S100a6, Rcn3, Serpinh1, Aebp1, Vim, Mmp2, Lsp1, Anxa1, Anxa2
## 4 Ccl5, Hp, H2-Ab1, Lyz2, Ccl6, H2-Aa, Ctss, Cd14, Atf3, Cd68, H2-Eb1, AW112010, S100a4, Cd24a, Fcer1g, Cfp, Trf, Aif1, Tyrobp, Ms4a6c, Tspan13, Cdkn1a, Ctsc, Tgm2, Grn
## 5 Ccl5, H2-Ab1, Lyz2, H2-Aa, Ctss, C1qc, Ctsd, Cd14, Cd68, H2-Eb1, AW112010, Fcer1g, Trf, Aif1, Tyrobp, Ms4a6c, Tspan13, Cdkn1a, Ctsc, Tgm2, Ctsb, Lgmn, Grn
## 6 Cxcl14, Apod, Mt1, Lcn2, Lpl, Clec3b, Gsn, Dcn, Serpinb1a, Igfbp4, Efemp1, Procr, Igfbp6, Col4a1, Rarres2, Col1a1, Col1a2, Fn1, Htra3, Ogn, Dpt, Ly6c1, Col3a1, Id3, Lum, Gadd45g, Ctsl, Col6a1, Mfap5, Sparc, Serping1, Glul, Col6a2, Gpx3, Ly6a, Dpep1, Bgn, Serpinf1, Col5a2, Ctsk, Pcolce, Mfap2, Loxl1, Fstl1, Phlda1, Plpp3, S100a6, Rcn3, Serpinh1, Aebp1, Pdpn, Dbi, C1s1, Mmp2, Anxa1, Anxa2
## 7 Actg2, Gng11, Acta2, Tagln, Mylk, Myl9, Cnn1, Col4a2, Postn, Tpm2, Tpm1, Col4a1, Rarres2, Col6a3, Col1a1, Col1a2, Fn1, S100a4, Col3a1, Crip2, Sparc, Csrp1, Mfge8, Bgn, Lgals1, Col5a2, Pcolce, Fstl1, Igfbp7, Cystm1, S100a6, Rcn3, Serpinh1, Aebp1, Vim, Anxa1, Anxa2
## 8 Cxcl14, Apod, Mt1, Dcn, Postn, Igfbp4, Igfbp6, Tpm1, Itm2a, Rarres2, Col6a3, Col1a1, Col1a2, Fn1, Ogn, Dpt, Col3a1, Rbp1, Lum, Gadd45g, Ctsl, Col6a1, Mfap5, Sparc, Serping1, Glul, Mfge8, Col6a2, Gpx3, Dpep1, Bgn, Lgals1, Tsc22d1, Serpinf1, Col5a2, Ctsk, Cnn3, Pcolce, Mfap2, Loxl1, Fstl1, Mmp14, Phlda1, S100a6, Rcn3, Serpinh1, Aebp1, Pdpn, Dbi, Vim, C1s1, Tmem176a, Mmp2, Lsp1, Anxa1, Anxa2
## 9 Ccl5, Nkg7, AW112010, Tspan13
## 10 H2-Ab1, Lyz2, Ccl6, H2-Aa, Pf4, Mt1, Ctss, C1qc, Ctsd, Cd14, Atf3, Cd68, H2-Eb1, AW112010, Fcer1g, Cfp, Trf, Aif1, Tyrobp, Ms4a6c, Fxyd2, Cdkn1a, Ctsc, Ctsb, Cystm1, Dbi, Lgmn, Tmem176a, Grn
## 11 H2-Ab1, Lyz2, H2-Aa, Ctss, C1qc, Ctsd, Cd14, Cd68, H2-Eb1, AW112010, Fcer1g, Cfp, Trf, Tuba1b, Aif1, Tyrobp, Ms4a6c, Cdkn1a, Ctsc, Tgm2, Ctsb, Lgmn, Grn
## 12 Actg2, Gng11, Acta2, Tagln, Myl9, Dcn, Postn, Tpm2, Tpm1, Col4a1, Col6a3, Col1a1, Col1a2, Fn1, S100a4, Col3a1, Tuba1b, Lum, Crip2, Col6a1, Sparc, Csrp1, Mfge8, Col6a2, Bgn, Lgals1, Serpinf1, Col5a2, Cnn3, Pcolce, Loxl1, Fstl1, Igfbp7, Phlda1, S100a6, Rcn3, Serpinh1, Vim, Anxa1, Anxa2
## cell_type celltype_score
## 1 Luminal Progenitor Cell 0.83
## 2 Luminal Progenitor Cell 0.88
## 3 Myoepithelial Cell 0.79
## 4 Luminal Progenitor Cell 0.79
## 5 Luminal Progenitor Cell 0.73
## 6 Luminal Progenitor Cell 0.81
## 7 Myoepithelial Cell 0.79
## 8 Luminal Progenitor Cell 0.76
## 9 Killer Cell, Luminal Cell 0.58, 0.58
## 10 Luminal Progenitor Cell 0.73
## 11 Luminal Progenitor Cell 0.73
## 12 Myoepithelial Cell 0.78
## celltype_related_marker
## 1 Aldh1a3, Plet1, Igfbp5, Glul, Anxa2, Phlda1
## 2 Kit, Lrg1, Slc12a2, Ptn, Tsc22d1, Anxa2
## 3 Acta2, Tagln, Actg2, Igfbp2, Cnn1, Myl9, Tpm2, Tpm1, Csrp1, Rbp1, Mfge8, Col4a2, Serpinh1
## 4 Cd14, Tgm2
## 5 Cd14, Tgm2
## 6 Glul, Lcn2, Dbi, Anxa2, Phlda1
## 7 Acta2, Tagln, Actg2, Cnn1, Mylk, Myl9, Tpm2, Tpm1, Csrp1, Mfge8, Col4a1, Col4a2, Serpinh1
## 8 Mfge8, Glul, Tsc22d1, Dbi, Anxa2, Phlda1
## 9 AW112010, Ccl5, Tspan13, Nkg7
## 10 Cd14, Dbi
## 11 Cd14, Tgm2
## 12 Acta2, Tagln, Actg2, Myl9, Tpm2, Tpm1, Csrp1, Mfge8, Col4a1, Serpinh1
## PMID
## 1 29225342, 29775597
## 2 21132007, 29225342, 29775597
## 3 29225342, 29775597
## 4 29225342, 29775597
## 5 29225342, 29775597
## 6 29775597
## 7 29225342, 29775597
## 8 29775597
## 9 29775597
## 10 29225342, 29775597
## 11 29225342, 29775597
## 12 29225342, 29775597Method 2- SingleR
mRNAseq.se <- MouseRNAseqData()## snapshotDate(): 2020-10-27## see ?celldex and browseVignettes('celldex') for documentation## loading from cache## see ?celldex and browseVignettes('celldex') for documentation## loading from cachemRNAseq.se## class: SummarizedExperiment
## dim: 21214 358
## metadata(0):
## assays(1): logcounts
## rownames(21214): Xkr4 Rp1 ... LOC100039574 LOC100039753
## rowData names(0):
## colnames(358): ERR525589Aligned ERR525592Aligned ... SRR1044043Aligned
## SRR1044044Aligned
## colData names(3): label.main label.fine label.ontmm10_ann_2 <- SingleR(test = scRNA_combined_mm10[['integrated']]@data,
ref = mRNAseq.se,
assay.type.test=1,
labels = mRNAseq.se$label.main)
mm10_ann_2## DataFrame with 4017 rows and 5 columns
## scores first.labels
## <matrix> <character>
## AAACCTGCATCACCCT-1_1 0.2749191: 0.0953161:-0.00695552:... Fibroblasts
## AAACCTGCATTGGGCC-1_1 0.1399648: 0.0770076: 0.09786590:... Endothelial cells
## AAACCTGGTGGTTTCA-1_1 0.3655346: 0.2011414: 0.07613095:... Fibroblasts
## AAAGCAAAGGATCGCA-1_1 0.0561536:-0.0161654: 0.16629786:... Monocytes
## AAAGCAAAGTCATGCT-1_1 0.0779979: 0.0166078: 0.10401372:... Monocytes
## ... ... ...
## TTTGGTTGTGACGGTA-1_4 0.37394650:0.1605742:0.0708835:... Fibroblasts
## TTTGGTTGTTGTCTTT-1_4 0.25604361:0.1082203:0.1364289:... Fibroblasts
## TTTGGTTTCTTTAGGG-1_4 0.30681930:0.1277484:0.0688820:... Fibroblasts
## TTTGTCAAGAGTGAGA-1_4 0.00969858:0.0060175:0.3195598:... T cells
## TTTGTCAAGCTTATCG-1_4 0.35637833:0.1532604:0.0670055:... Fibroblasts
## tuning.scores labels pruned.labels
## <DataFrame> <character> <character>
## AAACCTGCATCACCCT-1_1 0.326946:0.274919 Fibroblasts Fibroblasts
## AAACCTGCATTGGGCC-1_1 0.218958:0.144200 Endothelial cells Endothelial cells
## AAACCTGGTGGTTTCA-1_1 0.432832:0.365535 Fibroblasts Fibroblasts
## AAAGCAAAGGATCGCA-1_1 0.366443:0.228990 Monocytes Monocytes
## AAAGCAAAGTCATGCT-1_1 0.284693:0.222767 Monocytes Monocytes
## ... ... ... ...
## TTTGGTTGTGACGGTA-1_4 0.385029:0.07043289 Fibroblasts Fibroblasts
## TTTGGTTGTTGTCTTT-1_4 0.409715:0.00152361 Fibroblasts Fibroblasts
## TTTGGTTTCTTTAGGG-1_4 0.361933:0.30681930 Fibroblasts Fibroblasts
## TTTGTCAAGAGTGAGA-1_4 0.353053:0.28400529 T cells T cells
## TTTGTCAAGCTTATCG-1_4 0.437798:0.35637833 Fibroblasts FibroblastsSession info
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Spain.1252 LC_CTYPE=Spanish_Spain.1252
## [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Spain.1252
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] celldex_1.0.0 org.Mm.eg.db_3.12.0
## [3] org.Hs.eg.db_3.12.0 AnnotationDbi_1.52.0
## [5] clusterProfiler_3.18.1 fgsea_1.16.0
## [7] tibble_3.1.6 msigdbr_7.5.1
## [9] SingleR_1.4.1 SummarizedExperiment_1.20.0
## [11] Biobase_2.50.0 GenomicRanges_1.42.0
## [13] GenomeInfoDb_1.26.7 IRanges_2.24.1
## [15] S4Vectors_0.28.1 BiocGenerics_0.36.1
## [17] MatrixGenerics_1.2.1 matrixStats_0.62.0
## [19] BiocManager_1.30.17 ggplot2_3.3.5
## [21] scCATCH_3.2 purrr_0.3.4
## [23] sp_1.4-7 SeuratObject_4.1.0
## [25] Seurat_4.1.1 stringr_1.4.0
## [27] dplyr_1.0.8 R.utils_2.11.0
## [29] R.oo_1.24.0 R.methodsS3_1.8.1
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 reticulate_1.25
## [3] tidyselect_1.1.2 RSQLite_2.2.12
## [5] htmlwidgets_1.5.4 grid_4.0.4
## [7] BiocParallel_1.24.1 Rtsne_0.16
## [9] scatterpie_0.1.7 munsell_0.5.0
## [11] codetools_0.2-18 ica_1.0-2
## [13] future_1.25.0 miniUI_0.1.1.1
## [15] withr_2.5.0 spatstat.random_2.2-0
## [17] colorspace_2.0-3 GOSemSim_2.16.1
## [19] progressr_0.10.0 highr_0.9
## [21] knitr_1.31 rstudioapi_0.13
## [23] ROCR_1.0-11 tensor_1.5
## [25] DOSE_3.16.0 listenv_0.8.0
## [27] labeling_0.4.2 GenomeInfoDbData_1.2.4
## [29] polyclip_1.10-0 farver_2.1.0
## [31] bit64_4.0.5 downloader_0.4
## [33] parallelly_1.31.1 vctrs_0.4.1
## [35] generics_0.1.2 xfun_0.21
## [37] BiocFileCache_1.14.0 R6_2.5.1
## [39] graphlayouts_0.8.0 rsvd_1.0.5
## [41] bitops_1.0-7 spatstat.utils_2.3-1
## [43] cachem_1.0.4 DelayedArray_0.16.3
## [45] assertthat_0.2.1 promises_1.2.0.1
## [47] scales_1.2.0 ggraph_2.0.5
## [49] enrichplot_1.10.2 rgeos_0.5-9
## [51] gtable_0.3.0 beachmat_2.6.4
## [53] globals_0.15.0 goftest_1.2-3
## [55] tidygraph_1.2.1 rlang_1.0.2
## [57] splines_4.0.4 lazyeval_0.2.2
## [59] spatstat.geom_2.4-0 yaml_2.2.1
## [61] reshape2_1.4.4 abind_1.4-5
## [63] httpuv_1.6.5 qvalue_2.22.0
## [65] tools_4.0.4 ellipsis_0.3.2
## [67] spatstat.core_2.4-4 jquerylib_0.1.4
## [69] RColorBrewer_1.1-3 ggridges_0.5.3
## [71] Rcpp_1.0.8.3 plyr_1.8.7
## [73] sparseMatrixStats_1.2.1 progress_1.2.2
## [75] zlibbioc_1.36.0 RCurl_1.98-1.6
## [77] prettyunits_1.1.1 rpart_4.1-15
## [79] deldir_1.0-6 viridis_0.6.2
## [81] pbapply_1.5-0 cowplot_1.1.1
## [83] zoo_1.8-10 ggrepel_0.9.1
## [85] cluster_2.1.0 magrittr_2.0.1
## [87] RSpectra_0.16-1 data.table_1.14.2
## [89] scattermore_0.8 DO.db_2.9
## [91] lmtest_0.9-40 RANN_2.6.1
## [93] fitdistrplus_1.1-8 hms_1.1.1
## [95] patchwork_1.1.1 mime_0.12
## [97] evaluate_0.15 xtable_1.8-4
## [99] gridExtra_2.3 compiler_4.0.4
## [101] shadowtext_0.0.8 KernSmooth_2.23-18
## [103] crayon_1.5.1 htmltools_0.5.2
## [105] ggfun_0.0.6 mgcv_1.8-33
## [107] later_1.3.0 tidyr_1.2.0
## [109] DBI_1.1.2 ExperimentHub_1.16.1
## [111] tweenr_1.0.2 dbplyr_2.1.1
## [113] rappdirs_0.3.3 MASS_7.3-56
## [115] babelgene_22.3 Matrix_1.4-1
## [117] cli_3.2.0 igraph_1.3.0
## [119] pkgconfig_2.0.3 rvcheck_0.2.1
## [121] plotly_4.10.0 spatstat.sparse_2.1-1
## [123] bslib_0.3.1 XVector_0.30.0
## [125] yulab.utils_0.0.4 digest_0.6.27
## [127] sctransform_0.3.3 RcppAnnoy_0.0.19
## [129] spatstat.data_2.2-0 rmarkdown_2.14
## [131] leiden_0.4.2 fastmatch_1.1-3
## [133] uwot_0.1.11 DelayedMatrixStats_1.12.3
## [135] curl_4.3 shiny_1.7.1
## [137] lifecycle_1.0.1 nlme_3.1-152
## [139] jsonlite_1.8.0 BiocNeighbors_1.8.2
## [141] limma_3.46.0 viridisLite_0.4.0
## [143] fansi_1.0.3 pillar_1.7.0
## [145] lattice_0.20-41 fastmap_1.1.0
## [147] httr_1.4.3 survival_3.3-1
## [149] GO.db_3.12.1 interactiveDisplayBase_1.28.0
## [151] glue_1.4.2 png_0.1-7
## [153] BiocVersion_3.12.0 bit_4.0.4
## [155] ggforce_0.3.3 stringi_1.7.6
## [157] sass_0.4.1 blob_1.2.3
## [159] AnnotationHub_2.22.1 BiocSingular_1.6.0
## [161] memoise_2.0.1 irlba_2.3.5
## [163] future.apply_1.9.0